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Abstract 

A  non-isothermal,  non-isobaric  water  and  thermal  management  model  with  phase  change  has  been  developed  to  simulate  the  mass  and  energy 
(sensible  heat,  latent  heat,  chemical  reaction  energy,  electrical  energy)  transfer  processes  inside  a  PEM  fuel  cell  unit  with  a  non-uniform  stack 
temperature  (hereafter  the  stack  temperature  has  the  same  meaning  of  the  cell  temperature  since  only  single  cell  is  considered).  Based  on  this  model, 
the  following  parameters  can  be  predicted  along  the  channels  on  both  cathode  and  anode  sides:  current  density,  output  voltage,  stack  temperature, 
stream  pressure  and  temperatures,  stream  velocity,  relative  humidity,  water  vapor  mole  fraction,  water  liquid  fraction,  density,  viscosity,  Reynolds 
number,  required  pumping  power,  and  so  on.  Also,  the  results  from  uniform  stack  temperature  model  and  non-uniform  stack  temperature  model 
are  compared  in  this  research. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Most  current  automobiles  are  driven  by  internal  combustion 
engines  which  consume  fossil  fuel  and  generate  air  pollution. 
With  the  increasing  public  concerns  of  environmental  protec¬ 
tion,  it  is  predictable  that  more  and  more  strict  regulations  will  be 
enforced  to  reduce  or  limit  the  emission  of  these  vehicles  in  the 
future.  For  example,  California’s  zero  emission  vehicle  (ZEV) 
mandate  [1]  requires  10%  of  the  vehicles  sold  by  the  automotive 
manufacturers  after  year  2004  to  be  ZEVs  [2].  Similarly,  Euro¬ 
pean  auto-companies  are  required  to  meet  their  voluntary  carbon 
dioxide  (CO2)  emission  limits  set  by  the  European  Union  [3]. 
According  to  the  Kyoto  Protocol,  the  international  community 
is  committed  to  cut  greenhouse  gas  emissions  step  by  step.  The 
CO2  emitted  by  automobiles  is  one  of  the  major  components 
of  greenhouse  effects.  World  Governments,  including  the  Cana¬ 
dian  Government,  have  already  invested  a  lot  in  exploring  new 
ways  to  replace  the  internal  combustion  engine  in  automobiles. 
Among  all  the  technical  proposals,  the  fuel  cell  is  one  of  the 
most  potential  and  feasible  solutions  to  achieve  this  goal.  The 
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benefits  of  using  fuel  cells  are  as  follows  [4] :  firstly,  fuel  cells 
consume  hydrogen  instead  of  the  exhaustible  fossil  fuel,  which 
eventually  protect  our  natural  resource  and  environment;  sec¬ 
ondly,  fuel  cells  emit  only  water,  therefore,  there  is  no  pollution 
at  all.  Among  all  the  currently  existing  fuel  cells,  the  proton 
exchange  membrane  (PEM)  fuel  cell  has  been  widely  consid¬ 
ered  as  one  of  the  most  promising  candidates  for  automobiles 
since  it  has  one  additional  advantage  over  many  other  fuel  cells; 
the  PEM  can  operate  at  room  temperature  for  quick  startup. 

In  recent  years,  many  computer  models  of  PEM  fuel  cell 
have  been  developed.  Fuller  and  Newman  [5]  developed  a  non- 
isothermal  model  by  including  the  material  balances  in  the 
channel,  the  concentration  and  temperature  gradients  along  the 
channel  as  well  as  across  the  membrane  surface.  Nguyen  and 
White  [6]  studied  the  variation  in  current  density,  water  trans¬ 
port,  and  flow  temperatures  along  the  channel.  They  also  mod¬ 
eled  the  effect  of  varying  anode  inlet  humidity.  Subsequently, 
an  advanced  model  was  developed  by  Yi  and  Nguyen  [7]  to 
compare  different  fuel  cell  designs  with  coflow  and  counterflow 
heat  exchangers.  In  this  analysis,  they  include  the  thermal  mass 
of  the  stack,  account  for  the  impact  of  the  pressure  difference 
between  anode  and  cathode  on  water  transport  in  the  cell.  How¬ 
ever,  a  clear  validation  of  results  in  both  of  the  analyses  was 
not  presented.  Mosdale  and  Srinivasan  [8]  gave  a  good  review 
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Nomenclature 

a  water  vapor  activity  in  stream 

Ac  the  heat  transfer  area  between  stack  and  flow  in  a 
control  volume  (cm2) 

Acooi  the  heat  transfer  area  between  stack  and  coolant 
(cm2) 

Across  the  cross-section  area  of  channel  (cm2) 

As  the  cross-section  area  of  stack  (cm2) 

Cmw  concentration  of  water  at  interface  of  the  mem¬ 
brane  (mol  cm-3) 

Cpj  heat  capacity  of  species  i  (J  mol-1  K-1) 

d  channel  height  (cm) 

dp  pressure  drop  (Pa) 

D  hydraulic  diameter  of  channel  (cm) 

Dw  diffusion  coefficient  of  water  (cm2  s-1) 

D°  parameter  used  in  expression  for  diffusion  coeffi¬ 
cient  of  water  (cm2  s-1) 
fix)  friction  factor 

F  Faraday  constant,  96,487  C  equiv.-1 

h  channel  width  (cm) 

I  current  (A) 

I(x)  current  density  (A  cm-2) 

7°  exchange  current  density  for  the  oxygen  reaction 

(A  cm-2) 

k  thermal  conductivity  of  the  stream  (W  m-  1  °  C- 1 ) 

kc  evaporation  and  condensation  rate  constant  (s-1) 
L  length  of  channel  (cm) 

M  equivalent  weight  of  a  dry  membrane  (g  mol  - 1 ) 
Mi  molecular  weight  of  species  i  (g  mol-1) 
nd  electro-osmotic  drag  coefficient  (number  of  water 

molecules  carried  by  per  proton) 

N  mole  number  of  species  in  the  stream  (mol  s-1) 

Ach  number  of  channel  (s) 

NE  number  of  electrons  (A-1  s-1) 

Pi  partial  pressure  of  species  i  (Pa) 

P  cell  total  pressure  (Pa) 

^pump  pumping  power  (W) 

q  energy  (Js-1) 

Q  volume  flowrate  (m3  s-1) 

Ru  the  ideal-gas  constant  (8.3144  J  mol-1  K-1) 

Re  Reynolds  number 

RH  relative  humidity 

S  entropy  ( J  mol - 1  K- 1 ) 

t  thickness  (cm) 

T  temperature  of  stream  (K) 

Fs  temperature  of  stack  (K) 

U  overall  heat  transfer  coefficient 
(J  s-1  cm-2  °C-1) 

V  flow  velocity  (ms-1) 

VCeii  cell  voltage  (V) 

v  direction  along  the  channel  length 

y  direction  normal  to  the  channel  length 


Greek  letters 

a  excess  coefficients  of  air  (oxygen) 

o' area  reaction  area  coefficient 

/3h2  mole  fraction  of  hydrogen  (100%) 

/3o2  mole  fraction  of  oxygen  in  air  (20.9%) 

0  water  content  in  stream 

q  overpotential  for  the  oxygen  reaction  (V) 

li  dynamic  viscosity  (N  s  m-2) 

p  density  (kg  m- 3 ) 

Pm, dry  density  of  a  dry  membrane  (g  cm-3) 
crm  membrane  conductivity  (S3-1  cm-1) 

Subscripts 

air  dry  air 

avg  average 

A  anode 

cell  the  unit  fuel  cell 

concentration  the  concentration  of  species  in  the  streams 

cool  cooling  system 

C  cathode 

drag  electro-osmotic  drag 

e  electron 

ele  electricity 

gas  species  except  water  vapor  and  water  liquid 

gen  generated 

heat  heat 

H2  hydrogen 

H2O  produced  water 

in  inlet  of  channel 

latent  latent  heat 

liquid  liquid  water  in  the  flow 

loss  energy  released  to  environment 

m  membrane 

MW  water  in  membrane 

N2  nitrogen 

oc  open  circuit 

O2  oxygen 

pressure  the  partial  pressure  in  the  streams 

pump  the  pump  supplied  the  reactants  to  fuel  cell 

room  environment 

rxn  reaction 

sat  saturation 

sen  sensible 

stack  the  stack  of  fuel  cell 

system  a  closed  system 

vapor  water  vapor  in  the  flow 

water  all  water  including  vapor  and  liquid  in  the  flow 
0  in  the  standard  state 

1  A  per  ampere 

#  cathode  or  anode 
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on  the  work  that  has  been  conducted  at  Texas  A&M  Univer¬ 
sity.  In  this  review,  various  models  are  compared  and  the  effect 
of  different  humidification  on  the  performance  of  the  fuel  cell 
is  discussed.  Unfortunately,  this  paper  does  not  give  the  exact 
details  of  their  modeling  approaches.  Amphlett  et  al.  [9,10]  have 
developed  PEM  fuel  cell  models  that  are  based  on  both  the  the¬ 
oretical  mechanistic  analysis  and  empirical  data.  They  perform 
an  empirical  treatment  of  the  membrane  and  explain  the  water 
transport  processes  in  the  membrane.  Nevertheless,  this  model 
cannot  predict  some  parameters,  such  as  anode  humidification, 
temperature,  pressure  and  so  on.  Researchers  such  as  Marr  and 
Li  [11],  Dannenberg  et  al.  [12],  Hertwig  et  al.  [13],  Ge  and  Yi 
[14]  and  Xue  et  al.  [15]  have  been  performing  fuel  cell  modeling 
for  many  years  and  have  made  very  impressive  progress.  Those 
models  emphasized  important  characteristics  of  the  membrane, 
electrodes,  as  well  as  a  detailed  description  of  the  water  content 
in  the  membrane.  In  order  to  simplify  the  process,  most  models 
are  isothermal  and  isobaric  and  with  a  uniform  stack  tempera¬ 
ture.  These  parameters,  e.g.,  stream  temperature  and  pressure, 
stack  temperature  are  crucial  for  optimizing  PEM  fuel  cells. 
However,  few  papers  in  the  available  literature  addressed  them 
in  detail.  Zhou  et  al.  [16]  from  the  same  group  as  the  present 
paper,  proposed  a  non-isothermal  and  non-isobaric  model  with 
phase  change  effects  to  address  these  challenges,  although  the 
stack  temperature  was  also  assumed  to  be  uniform  to  simplify 
the  real  complex  processes  inside  the  PEM  fuel  cells. 

In  the  present  study,  the  model,  proposed  by  Zhou  et  al.  [16], 
will  be  extended  to  a  non-isothermal  and  non-isobaric  model 
with  non-uniform  stack  temperature.  This  improved  model  and 
some  simulation  results  are  presented  below. 

2.  Mathematical  model 

Fig.  1  shows  a  typical  construction  of  a  PEM  fuel  cell.  In  the 
present  study,  two  coordinate  axes  are  chosen.  The  v- axis  is  along 
to  the  gas  channels.  The  temperature,  pressure  and  concentration 
of  gas  flow  will  be  calculated  along  this  direction.  The  y-axis  is 
perpendicular  to  the  membrane.  The  hydrogen  ions  and  water 
molecules  transport  from  anode  to  cathode  along  this  direction. 


Membrane 


Fig.  1.  Schematic  diagram  of  PEM  fuel  cell  modeling  regions. 


The  model  includes  the  following  parts:  mass  balance,  energy 
balance,  pressure  drop  and  cell  output  voltage. 

Although  most  of  the  equations  have  been  presented  in  Zhou 
et  al.  [16],  they  are  repeated  here  for  completeness. 

2.7.  Basic  assumptions 

In  the  present  model,  the  water  transport  mechanisms  in  the 
membrane  are  based  on  the  study  by  Yi  and  Nguyen’s  study  [7]. 
Some  corresponding  assumptions  are  listed  as  follows: 

1.  Only  water  vapor  can  diffuse  into  the  electrode  and  pass 
through  the  membrane. 

2.  The  electrode  layer  is  “ultra  thin”,  gas  diffusion  through  the 
electrode  porous  layer  is  neglected. 

3.  The  gases  and  water  vapor  are  fully  mixed,  and  the  mixture 
is  ideal  gas. 

4.  Liquid  water  exists  only  in  the  form  of  small  droplets  and  the 
volume  is  negligible. 

5 .  Water  vapor  is  produced  in  the  electrochemical  reaction.  This 
assumption  could  bring  some  inaccuracy  since  the  operating 
temperature  is  usually  under  100  °C.  In  reality,  the  product 
water  could  be  liquid  water,  or  part  of  water  vapor  and  part 
of  liquid  water  depending  on  the  local  conditions. 

6.  No  voltage  drop  exists  along  the  flow  channels.  This  assump¬ 
tion  was  made  to  simplify  the  real  situation.  It  is  a  reasonable 
assumption  because  the  bipolar  plate,  usually  made  from 
graphite,  has  a  high  electrical  conductivity. 

7.  The  channels  in  each  unit  cell  have  the  same  geometry  and 
same  surface  roughness. 

8.  A  single  channel  is  assumed  to  represent  the  unit  cell  for 
numerical  simulations. 

9.  The  temperature  of  the  solid  (including  ME  A  and  plates)  is 
assumed  to  be  constant  in  the  y-direction  only. 

2.2.  Mass  balance 

Fig.  2  shows  the  mass  balance  in  a  unit  fuel  cell.  The  amount 
of  inlet  gases  is  calculated  according  to  the  amount  of  gases 
consumed  by  the  electrochemical  reaction  for  PEM  fuel  cell: 
2H2  +  02  =  2H2O.  One  equivalence  of  electrons  is  1  mole  of 
electrons  or  6.022  x  1023  electrons  (Avagadro’s  number).  This 
quantity  of  electrons  has  the  charge  of  96,487  coulombs  (C) 
(Faraday’s  Constant).  Therefore,  the  charge  of  a  single  electron 

Air,vapor  H2,  vapor 
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is  1.602  x  10  19  C.  One  ampere  of  current  is  defined  as  1  C  s  1 . 
As  a  result,  the  mole  number  of  electrons  for  1  A  of  current  is 
as  follows: 

«e  = - NElA  „  =  1.03656546  x  10-5  mole  (As)-1  (1) 

6.022  x  1023 

where  NEi  a  represents  the  number  of  electrons  of  1  A.  It  can  be 
calculated  by  NEia  =  1/1.602  x  10~19.  Based  on  the  reaction 
equation,  the  theoretical  mole  numbers  of  consumed  oxygen 
and  hydrogen,  and  produced  water  for  1  A  current  output  can  be 
obtained  by  the  following  equations: 

«02,1A  =  \ne  (2) 

2  1 

^H20,1A  =  4  =  2ne  (3) 

^H2,1A  =  fwe  =  \nQ  (4) 

The  excess  coefficients  for  air  (oxygen)  or  hydrogen  is 
defined  as 


where  c^area  is  the  reaction  area  coefficient  that  accounts  for  the 
land  area  for  reaction  due  to  gas  diffusion  from  the  channel  to 
diffusion  layer. 

The  variations  of  vapor  and  liquid  water  along  the  chan¬ 
nels  are  more  complicated.  Water  vapor  transport  and  con¬ 
densation,  and  liquid  water  evaporation  are  considered  in  this 
model.  There  are  three  water  transport  mechanisms  across  the 
membrane,  according  Yi  and  Nguyen  [7]:  (a)  electro-osmotic 
drag — since  the  hydrogen  ions  pass  through  the  membrane,  the 
water  molecules  are  carried  from  the  anode  to  the  cathode;  (b) 
back-diffusion  by  the  concentration  gradient  of  water — because 
the  water  concentration  is  different,  some  water  molecules  dif¬ 
fuse  from  the  cathode  to  the  anode;  (c)  convection  by  the  pressure 
gradient — water  moves  from  higher-pressure  side  to  the  lower 
one.  In  a  calculated  volume,  usually  it  is  very  small  if  the  grid 
size  is  small  enough,  there  is  negligible  potential  gradient  in  the 
v-direction  within  the  calculated  volume.  The  electro-osmotic 
drag  flux  (mol  s 1)  in  the  y-direction  is  as  follows  [6]: 


a  = 


actually  supplied  mole  number  of  air  (oxygen)  or  hydrogen 
theoretically  consumed  mole  number  of  air  (oxygen)  or  hydrogen 


(5) 


Therefore,  the  supplied  oxygen,  nitrogen  and  hydrogen  mole 
numbers  can  be  calculated  by  the  following  equations: 


«02,in,l  A  =  ^02^02,l  A 


^N2,in,  1A  —  ^02,in,lA 


1  Po2 
Po2 


«H2,in,l  A 


^H2^H2,1  A 

Pu2 


(6) 

(7) 

(8) 


/^drag(};)  — 


nd(x)I(x) 

F 


(15) 


where  I{x)  is  the  local  current  density  of  the  fuel  cell,  F  is  Fara¬ 
day’s  constant,  nd  the  electro-osmotic  drag  coefficient  which 
represents  the  number  of  water  molecules  carried  by  one  proton 
and  is  calculated  by  [6] : 


nd(x)  = 


0.0049  +  2.02 a(x)  —  4.53 a2(x)  +  4.09 a3(x)  0  <  a(x)  <  1 
1.5849  +  0.159(a(jc)  -  1)  a(x)  >  1 


(16) 


where  po2  Is  the  mole  fraction  of  oxygen  in  air  (Po2  =  20.9%) 
and  /3h2  is  that  of  the  hydrogen  in  anode  (/3h2  =  1  in  the  present 
study). 

For  generating  I  amperes  of  current,  the  molar  flowrates  of 
oxygen,  nitrogen  and  hydrogen  for  the  single  channel  are  eval¬ 
uated  in  (mol  s_1)  as 


/Vc,02,in  = 


7^Q2,in,l  A 


Nch 


/Vc,N2,in  =  /Vc,02,in 


i  -  A> 


PO: 


NA,R2,in  =  I 


«H2,in,l  A 

Nch 


(9) 

(10) 

(11) 


The  components  of  the  mixture  vary  along  the  channels  and 
the  local  molar  flowrates  in  channel  are  defined  as  follows: 


d/Vc,o2(*) 

dx 

&Nc,n2(x) 

dx 

dAA,H2M 

dx 


—  no2,\Al(x)haaiea 


=  0 


=  -nu2,iAI(x)ha 


area 


(12) 

(13) 

(14) 


where  a(x)  is  the  activity  of  water  vapor  in  membrane.  Since 
the  membrane  is  placed  between  the  cathode  and  anode,  the 
water  vapor  activity  in  the  membrane  is  affected  by  the  water 
vapor  activity  at  both  cathode  and  anode.  A  weighted  average 
water  activity  for  the  water  vapor  activity  in  the  membrane  is 
employed.  The  water  vapor  activity  at  anode  or  cathode  a#(x )  is 
defined  as  [6] : 

,  x  /V#, vapor C*0  P#(x) 

a#(x)  =  — — - - 

ZiN#Ax)  P#, satto 

(/,  the  species  in  the  flow  stream#)  (17) 

The  water  vapor  activity  at  the  membrane  is  defined  as 
aM(x)  =  aMaA(x)  +  (1  -  aM)ac(x)  (18) 


where  c^m  is  the  weight  coefficient.  Thus,  the  water  vapor  activity 
at  the  membrane  depends  on  the  water  activity  at  both  cathode 
and  anode. 

The  diffusion  flux  caused  by  the  concentration  gradient  of 
water  can  be  written  as  follows: 


/^concentration  (}0  —  ^MW 


(  9cmw 


(19) 
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where  the  diffusion  coefficient  of  water  is  as  follows  [7] : 


£>mw  =  < 


(0.0049  +  2.02 a(x)  -  4.53a2(x)  +  4.09a3 (x))D°  exp 

1  1 


2416 


1 


1 


303  Ts(x) 


1.59  +  0.159 [(a(x)~  l)]D°exp 


2416 


303  Ts(x ) 


for  0  <  a(x)  <  1 


for  a(x)  >  1 


(20) 


Convection  flux  caused  by  pressure  gradient  [7]  is 


pressure 


(y)  = 


CMW.cW  +  CmW.aW 


kp  f  dp# 

, vapor  \ 

) 


MwW 


(21) 


where  p#, vapor  is  the  water  vapor  pressure  in  the  anode  and  cathode  channels,  respectively.  kp  the  permeability  of  water  in  the 
membrane,  gw  M  the  water  viscosity  and  cmw,cW  and  cmw,aM  are  the  concentrations  of  water  in  cathode  and  anode,  respectively, 
with  the  expression  given  below  [6]: 

f  Pm,  dry 


CMW,#W  =  < 


37m,  dry 
Pm,  dry 


[0.043  +  17.8a#(v)  —  39.8 a#(x)  +  36.0 a#(x)]  for  0  <  a#(x)  <  1 


M 


[14  +  1.4(a#(x)  —  1)] 


for  a#(x)  >  1 


(22) 


m,dry 


where  pm,dry  and  Mm^ ry  are  the  density  and  the  equivalent  weight  of  a  dry  proton  exchange  membrane.  Therefore,  the  change  of 
water  flux  along  the  channel  at  the  cathode  can  be  expressed  by 


fl^C,  water  CO 

dx 


—  [^H2o7(v)  +  N drag(}7)  ^concentration  (y)  ^pressure  (}7)]^^area 


(23) 


The  variable  trend  of  water  flux  in  anode  can  be  expressed  by 

water  00 


dx 


—  [  3/ drag(.y)  +  3/ concentration  (+)  +  N pressure  (}7)]^<^area 


(24) 


The  mole  number  of  water  condensation  or  evaporation  can  be  calculated  as  [6] : 


fl^#,  liquid  (a) 

dx 


krhd  \ 


.  ^#,vaporOO 
RuT#(x)J  L £/##,/(*) 


p#(x)  -  p#f satW 


(/,  the  gaseous  species  in  stream  #) 


(25) 


In  order  to  present  the  state  of  vapor  water  and  liquid  water,  relative  humidity  (RH)  and  relative  water  content  (0)  are  defined  as 
follows: 


partial  pressure  of  vapor  water  N#fWspor(x)  p#(x) 

Kri#(V)  =  - - — -  or  Kri#(V)  = 


saturation  pressure 


y ^jN#mj(x)  p#,  sat  (a) 


(26) 


and 


mole  number  of  water  (vapor  +  liquid)  V#,  waterO)  P#(x)  ...  .  . 

0#(jc)  =  - ; - ; - - - : — : - : - - -  or  0#(v)  =  — — - —  (i,  the  gaseous  species  in  stream#)  (27) 


mole  number  of  water  in  saturation 


ZiN#Ax)  P#, satto 


2.3.  Energy  balance 

As  shown  in  Fig.  3,  the  energy  in  a  unit  fuel  cell  consists  of  the  energy  released  from  chemical  reactions,  which  is  the  source  of 
energy  in  fuel  cell;  the  electrical  energy  for  generating  power;  the  heat  for  increasing  the  solid  phase  temperature;  the  sensible  heat 
for  increasing  the  temperature  of  flow;  the  latent  heat  for  water  phase  changes;  the  waste  heat  taken  away  by  coolant;  and  the  heat 
loss  to  the  environment  at  the  inlet  and  exit  of  channel. 

The  total  energy  balance  is 

%xn(x)  —  ^elecOO  T  ^heatOO  —  ^elecOO  T  #C*0stack  3“  ^  ^  ff#,sen(-K)  +  ^  ^  ff#, latent (x)  +  ^cool(A)  (28) 

The  energy  released  from  electrochemical  reaction  is  difficult  to  calculate.  In  general,  enthalpy  or  entropy  is  employed  to  evaluate 
the  energy  and  electrical  work  in  an  electrochemical  system.  For  a  reversible  cell: 

4heat«  =  T0AS(x)  (29) 

In  this  study,  the  entropy  is  used  to  calculate  the  released  energy.  For  the  reacting  and  non-reacting  system,  the  entropy  balance 
for  undergoing  any  process  can  be  expressed  as 

(5in  —  Sout)  +  S gen  —  ^4  ■Ssystem 


(30) 
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Sensible  heat  at  cathode  Sensible  heat  at  anode 


Fig.  3.  Energy  balance  of  a  unit  fuel  cell. 


released  heat.  The  total  reaction  heat  flux  can  be  expressed  as 


^heat(x)  — 


ASa(*)  AScW 


2  F 


4  F 


^s^(x)/z<Xarea  dx 


-  ^(x)/(x)/zaarea  dx 


(35) 


If  the  numerical  value  of  gheat(x)  is  negative,  it  means  that  the 
chemical  reaction  emits  heat  to  the  surrounding. 

The  electrical  energy  (J  s-1)  is 

tfelecCO  =  Kell(x)/(x)/lQ'areadx  (36) 


This  means  that  the  change  of  entropy  in  a  system  can  be  deter¬ 
mined  by  the  net  entropy  transfer  and  the  entropy  generated  in 
the  system.  The  difference  of  the  entropy  change  between  a  sys¬ 
tem  with  chemical  reaction  and  a  non-reacting  system  is  that: 
the  entropy  relations  for  the  reactants  and  the  products  involve 
the  entropies  of  the  components,  not  entropy  changes,  which  is 
the  case  for  non-reacting  system  [17].  Thus,  a  common  base  for 
the  entropy  of  all  substances  is  established  by  the  third  law  of 
thermodynamics.  Based  on  the  third  law  of  thermodynamics,  the 
entropy  has  a  common  universal  scale  for  each  chemical  com¬ 
pound.  The  common  scale,  called  the  absolute  entropy,  is  based 
on  the  fact  that  the  entropy  for  any  pure  element  is  zero  at  the 
temperature  of  absolute  zero.  For  the  electrochemical  reaction 
in  a  PEM  system,  it  is  difficult  to  calculate  the  reaction  heat  for 
the  total  reaction.  However,  if  one  could  tell  the  heat  generation 
in  which  electrode  heat  is  generated  and  in  which  it  is  absorbed, 
the  corresponding  total  reaction  heat  can  be  calculated  easily. 
For  the  electrode  reaction  at  the  anode,  the  entropy  change  can 
be  calculated  by 

=  4^h+  +  4Sc~  -  2Su2  (31) 

where  S°  is  the  absolute  values  of  species  at  the  standard  state: 
To  =  298.15  K,  Po  =  1  bar.  The  numerical  values  for  the  species 
are  taken  from  Refs.  [18,19].  The  entropy  change  at  anode  is: 
Sa  =  0.208  J  (mol  K)“ 1  [  1 8] .  In  the  same  way,  the  entropy  change 
for  the  cathode  reaction,  can  be  calculated  by 


A  Sc  =  2S§2o  -  S£2  -  4S°_  -  4S°  +  (32) 

The  entropy  change  at  the  cathode  is: 
Sc  =  —326.36  J  (mol  K)-1  [18].  The  absolute  entropy  of 
species  i  at  temperature  T  and  pressure  P  can  be  calculated 
from 

Si(T,  P,  x)  =  S°  +  [  [ 

4  To  *  JP0 

(33) 


dvi(x) 

dT(x) 


where  v/(x)  is  the  specific  volume  of  flow  at  the  location  of  x. 
Considering  the  total  reaction  in  a  cell,  the  reaction  heat  flux 
(J  s_1)  can  be  calculated  for  the  reversible  process  as 


^heat(x)  — 


ASa(x)  ASc(x) 
IF  4  F 


Ts(x)I(x)ha3irQaL  dx 


(34) 


For  the  real  irreversible  process,  due  to  the  ohmic  loss  and 
reaction  resistance,  some  of  the  electrical  energy  is  turn  into  the 


The  sensible  heat  of  mixture  flow  at  cathode  channel  (J  s  1 )  is 

dc[#,  sen  W  =  Y^^N#3x)Cpj(x)]dT#(x) 

i 

(/,  the  gaseous  species  in  stream#)  (37) 


The  latent  heat  for  water  vapor  condensation  or  liquid  water 
evaporation  at  cathode  or  anode  channels  (J  s-1): 

dc[#, latent (x)  =  [7/#?Va por(x)  H# , liqUid (x)]  d7V#Jiquid(x)  (38) 

In  this  model,  the  heat  taken  by  the  coolant  is  considered  and 
expressed  by 


^cool(-X)  —  U cod  AC001  [  7g  (a)  Tcooi] 


(39) 


where  Ucoo\  is  the  heat  transfer  coefficient  between  stack  and 
coolant,  Acooi  is  the  area  of  heat  transfer. 

When  the  streams  flow  along  the  channels,  they  will  gain  or 
lose  heat  due  to  the  heat  transfer  between  the  fluid  and  stack. 
Therefore,  the  temperature  of  flow  depends  on  the  stack  tem¬ 
perature  and  the  latent  heat  as  well,  which  can  be  calculated  as 
follows: 


53  [A t#Ax)cpAx)] 


moo 

dx 


—  [TTw,  vapor  (x) 


#W,liquid(x)] 


d  A/#  ?  iiqUid  (x) 
dx 


+  U#h[Ts(x)  -  7#(x)] 


(40) 


U #  represents  the  heat  transfer  coefficient  between  the  flow 
stream  #  and  stack.  The  term  on  left  side  of  the  equation  repre¬ 
sents  the  heat  flux  obtained  by  gaseous  flows.  The  first  term  on 
the  right  side  of  the  equation  accounts  for  the  enthalpy  change 
due  to  condensation  or  evaporation  of  water  in  the  channels, 
which  can  be  calculated  using  Eq.  (41).  The  second  term  on  the 
right  side  of  the  equation  is  for  the  convection  heat  flux  between 
the  stream  and  the  stack: 


H#,  vapor  (x)  FI#  liquid  (x) 

=  45070  -  41.9[7#(x)  -  273]  +  3.44  x  10"3[7#(x)  -  273]2 
+  2.54  x  10_6[7#(jc)  -  273]3 
-  8.98  x  10"10[7#(V>  -  273]4 


(41) 
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The  stack  temperature  varies  with  position/location  along  the 
channels.  It  is  determined  by  the  local  current  density,  local  latent 
heat,  and  local  cooling  heat.  It  is  calculated  by  the  following 
energy  balance  equation: 


d27s(x)  dTs(x) 

Ask  ,  ,  I  Vqijqnjji  A')  V'Vliquid  (-V )  |  Cp  ,  water  ^ 

=  -AcU#[Ta(x)  +  Tc(x)  -  27s  W] 


^cool  U cool  [Tcool  00  7s  W] 


[77a,  vapor  00  7/a, liquid^)] 


[77c,  vapor  00  7/c,liquidOO] 

(  A^aW  ,  AStOOO 


dAOyiiquidOO 

dx 

dAfc,  liquid  00 


dx 


+ 


V 


IF 


+ 


4  F 


■ 


Ts(x)  -  rj(x) 


I(x)ha 


area 


(42) 


where  k  is  the  thermal  conductivity  of  the  stack,  As  the  cross- 
section  area  of  the  stack  along  the  flow  direction,  Ac  the  heat 
transfer  area  between  stack  and  flows,  and  Acooi  is  the  heat  trans¬ 
fer  area  between  stack  and  coolant.  The  term  on  the  left  side  of 
equation  represents  energy  flow  by  the  conduction  in  the  stack  of 
the  cell  along  the  gas  flow  path  (v-direction).  The  temperature 
distribution  normal  to  gas  flow  (y-direction)  is  assumed  to  be 
uniform.  The  first  term  on  the  right  side  of  equation  is  for  con¬ 
vective  heat  transfer  between  the  streams  in  the  channels  and  the 
cell  stack.  The  second  term  on  the  right  side  of  equation  accounts 
for  convective  heat  transfer  between  the  cell  stack  and  coolant. 
The  temperature  of  coolant  is  assumed  to  be  constant  along  the 
channels.  The  third  and  fourth  terms  represent  the  energy  taken 
or  released  from  the  phase  changes  of  water  in  the  anode  or  cath¬ 
ode  flow  (latent  heat),  which  can  be  calculated  in  Eq.  (40).  The 
last  term  represents  heat  generation  by  the  reversible  chemical 
reaction  process. 


2.4.  Pressure  drop 

The  pressure  drop  of  the  gas  mixture  in  the  fuel  cell  flow 
channel  was  rarely  considered  in  currently  available  fuel  cell 
research  publications.  But,  in  industrial  design  and  practice,  it  is 
a  significant  parameter  simply  because  it  directly  affects  system 
efficiency. 

The  saturation  pressure  (Pa)  can  be  expressed  in  terms  of  the 
local  temperature  [19]: 


The  local  velocity  (ms  1 )  in  the  cathode  and  anode  can  be 
calculated  as  follows: 


V#(x)  = 


2#  00 


, cross 


(45) 


where  A#?cross  is  the  cross-section  area  of  channel. 

Since  the  mole  fraction  of  gases  in  the  channel  varies,  local 
density  (kgm-3)  also  varies  with  the  different  components  in 
the  flow.  It  can  be  calculated  from 


P#(x)  = 

i 


N#,i(x) 


P#00 

T#(x)Ru 


(/,  the  gaseous  species  in  stream  #) 


(46) 


Local  dynamic  viscosity  can  be  calculated  by  interpola¬ 
tion.  /x/,ioo  is  the  gas  dynamic  viscosity  at  100  °C  and  /x*,o  is 
the  gas  dynamic  viscosity  at  0°C.  The  temperature  range  of 
flow  in  the  calculated  cases  is  0-100  °C,  so,  the  local  dynamic 
viscosity  is 


m(x)  =  E  { 


AM*)  T#(x)  -  273 

E ,-AM*)  i  ioo  -  o 


(Mi.  ioo~Mi.o)  +  Mi.o 


(/,  the  gaseous  species  in  stream  #) 


(47) 


For  laminar  flow,  pressure  drop  in  each  control  volume  can 
be  expressed  as  (Pa): 

Vl(x) 


d  P#(x) 
dx 


=  P#(x)f#(x) 


J _ 

2D 


(48) 


where  f#(x)  is  the  friction  factor  and  D  the  hydraulic  diameter  of 
the  channel.  In  this  model,  the  channels  are  straight,  so  that  only 
friction  loss  is  considered.  The  local  pressure  (Pa)  is  calculated 
from  the  pressure  at  inlet  by  subtracting  the  pressure  drop  from 
the  inlet  to  the  current  control  volume: 


P#(X )  =  p#t in  - 


r0 


dp#  (x) 

dx 


dx 


(49) 


The  total  required  pumping  power  (W)  is  used  by  the  designer 
to  choose  a  pump  to  maintain  the  flow.  It  is  given  by 


*L 


ft, 


pump 


—  A^ch 


d /?#(*) 

dx 


Q#(x)  dx 


(50) 


Psai(*)=  1-013  x  105  x  io-2'1794+0-02953(7#W-273)-9-1837x10“5(7#W-173)2+1-4454x10“7(7#W-273)3 


(43) 


Based  on  the  assumption  that  the  mixture  is  regarded  as  an 
ideal  gas,  local  volumetric  flow  (m3s_1)  in  the  cathode  and 
anode  can  be  calculated  using  the  ideal-gas  law: 


Q#(x)  =  yN#j(x)Ru 

i 


T#(x) 

P#(x) 


(/,  the  gaseous  species  in  stream  #) 


(44) 


2.5.  Cell  output  voltage 

Cell  output  voltage  is  a  significant  parameter  that  is  used 
to  evaluate  the  performance  of  a  PEM  fuel  cell.  The  goal  in 
this  section  is  to  model  the  potential  losses  in  the  gas  diffusion 
layers  and  membrane,  so  that  the  output  potential  can  be  accu¬ 
rately  predicted.  The  output  voltage  of  the  fuel  cell  is  modeled  as 
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1  2  3 


k 


n-1  1 


inlet 

outlet 

Fig.  4.  Control  volume  layout  for  numerical  solution  along  the  channel. 


the  reversible  cell  voltage  minus  activation  losses,  concentration 
over-potential  and  ohmic  resistance  of  the  electrodes,  catalyst 
layer  and  membrane.  The  cell  potential  is  expressed  as  [6]: 

Well  =  v0c  -  rj(x)  -  (51) 

where  V0c  is  the  open-circuit  potential  of  fuel  cell  and  77 (v)  refers 
to  the  cell  over-potential  which  combines  activation  losses  and 
concentration  losses  together  and  is  calculated  as  follows  [6]: 


RuTs(x) 

0.5  F 


In 


m  \ 

I°Pc,o2(x)J 


(52) 


where  7°  is  the  exchange  current  density  at  one  atmosphere 
of  oxygen,  pc,o2(x)  the  partial  pressure  of  oxygen  at  cathode, 

1{x]tn)  the  ohmic  loss  which  depends  on  the  membrane  water 

am  OO  r 

activity,  stack  temperature  and  membrane  thickness  and  am(x) 
is  the  membrane  conductivity  calculated  by  the  following 
equation  [19]: 


<ym(x)  =  (  0.005 14  Mm,dry  cMw(>)  -  0.00326  ) 

Pm,  dry  / 


xexp  1268 


1 


1 


303  Ts(x ) 


(53) 


where 


Pm’dIL [0.043  +  17.8a(x)  -  39.85a2(x)  +  36.0a3  (x)] 


cmwM  =  < 


3Tm,dry 

Pm,  dry 


M 


[14+  \A(a(x)-  1)] 


m,dry 


The  current  density  is  distributed  unevenly  along  the  chan¬ 
nel.  However,  when  a  fuel  cell  is  evaluated,  the  output  current 
is  of  major  concern.  Therefore,  the  average  current  density  is 
calculated  using: 

i  p 

4vg  =  y  I(x)dx  (55) 

L  Jo 


3.  Solution  procedure 


Based  on  previous  analysis,  the  model  can  be  summarized 
into  a  set  of  non-linear  differential  equations.  In  this  study,  iter¬ 
ative  methods  are  chosen  to  solve  the  differential  equations.  The 
solution  procedure  includes  three  loops.  The  outer  loop  adopts 
a  Gauss-Jacobi  method  to  solve  the  equations,  such  as  mass 
of  hydrogen,  mass  of  oxygen,  pressure  drop,  stream  tempera¬ 
tures,  mass  of  water  vapor  and  liquid  water,  and  so  on.  In  the 
inner  loop,  a  Tri-Diagonal  Matrix  Algorithm  is  used  to  calculate 
the  energy  balance  to  get  the  stack  temperature.  For  the  non¬ 
linear  algorithm  equation,  such  as  cell  potential-current  density 


Table  1 


Geometry  parameters  of  a  single  fuel  cell. 


Parameters 

Values 

Channel  length  (L) 

85  cm 

Channel  width  at  cathode  and  anode  (h) 

0.15  cm 

Channel  height  at  cathode  and  anode  (< d) 

0.08  cm 

Channel  number  of  cathode  and  anode  (Nc\f) 

6 

The  effective  area 

100  cm2 

Condensation  rate  constant  ( kc ) 

1.0s-1 

Membrane  dry  density  (pm,dry) 

2.0  g  cm-3 

Membrane  dry  equivalent  weight  (Mm^ry) 

1 100  g  mol-1 

Membrane  thickness  (tm) 

0.01275  cm 

Fuel  cell  open-circuit  voltage  (Vqc) 

1.1  V 

Oxygen  exchange  current  density  (1°) 

0.01  A  cm-2 

Diffusion  coefficient  of  water  in  membrane  ( D° ) 

5.5  x  10-7  cm2  s-1 

Eq.  (51),  the  Newton-Raphson  method  is 

applied  to  calculate 

the  local  current  density.  Before  the  iterative  methods  can  be 
applied,  the  differential  equations  are  required  to  be  discretized 
into  algebraic  equations.  In  this  simulation,  the  finite  difference 
method  is  adopted.  The  channel  is  subdivided  into  n  control 
volumes  of  equal  length  L  =  Lin  in  v-direction  (shown  in  Fig.  4). 
The  exit  values  at  the  kth  control  volume  are  the  inlet  values  at 
the  ( k+  l)th  control  volume,  and  all  variables  are  stored  at  the 

for  0  <  a(x)  <  1 

(54) 

for  a(x)  >  1 

centroid  of  each  cell.  A  set  of  differential  equations  is  replaced 
by  algebraic  equations  based  on  this  method. 

In  this  simulation,  the  hydrogen  and  air  flow  in  the  channels 
are  in  the  coflow  mode.  Table  1  lists  the  basic  geometric  parame¬ 
ters  and  electrode  and  membrane  properties  of  the  unit  fuel  cell. 

4.  Model  validation  and  results  analysis 

A  Nexa™  system  [20]  is  used  in  this  experiment.  The 
Nexa™  system  consumes  hydrogen  and  air  to  provide  dc  power 
up  to  1200  W  with  a  nominal  output  voltage  of  26  VDC.  It  con¬ 
tains  a  BALLARD®  fuel  cell  stack,  as  well  as  all  the  auxiliary 
equipment  necessary  for  fuel  cell  operation. 

Fig.  5  compares  the  simulation  results  for  constant  stack 
temperature  and  for  variable  stack  temperature  with  the  experi¬ 
mental  data.  When  the  current  density  is  less  than  0.25  A  cm-2, 
the  output  voltage  in  the  modes  is  greater  than  the  experimental 
results.  After  the  current  density  exceeds  0.3  A  cm-2,  the 
output  voltage  in  the  model  is  less  than  the  experiment  results. 
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Fig.  5.  Comparison  of  the  model  predictions  with  the  empirical  data. 


Obviously,  the  inner  resistances  in  the  models  are  greater  than 
that  in  practice.  Meanwhile,  the  resistance  in  the  constant  stack 
temperature  model  is  greater  than  that  in  the  variable  stack 
temperature  model.  It  is  difficult  to  explain  clearly  the  exact 
reasons  for  the  errors  between  the  model  results  and  experiment 
results.  However,  some  factors  may  contribute  to  the  error  of  the 
resistance  calculation,  which  lead  to  the  above  simulation  result. 

In  the  model,  represents  the  ohmic  loss  inside  the 

fuel  cell,  which  depends  on  local  current  density,  membrane 
thickness,  and  the  membrane  conductivity.  The  experimental 
equation  for  membrane  conductivity  comes  from  experimental 
results  of  Ref.  [21],  which  is  based  on  a  fully  hydrated  mem¬ 
brane  of  Nation  1 17.  In  the  real  situation,  at  many  local  points  of 
membrane,  the  hydration  is  not  perfect,  which  means  the  water 
activity  in  membrane  varies.  In  addition,  at  the  time  of  this  report, 
data  on  the  material  and  membrane  thickness  for  the  Nexa™ 
power  module  was  not  available.  According  to  previous  research 
[22] ,  the  material  and  thickness  of  the  membrane  affect  the  mem¬ 
brane  conductivity  greatly.  Without  the  detailed  information 
regarding  the  membrane,  it  is  challenging  to  obtain  agreeable 
curves  between  model  results  and  experimental  results. 

4.1.  Base  case  analysis 

The  simulation  is  based  on  an  operating  condition  at  near 
room  temperature  and  low  pressure,  which  is  called  the  “Base 
case”.  The  detailed  operating  conditions  of  the  base  case  are 
listed  in  Table  2. 

Fig.  6  shows  the  local  current  density  in  the  base  case  along 
the  channel.  The  current  density  is  highest  at  the  channel  inlet.  It 
then  drops  quickly  to  the  lowest  point.  After  that,  current  density 
increases  slowly  along  the  channel  until  it  reaches  the  channel 
exit,  where  the  current  density  increases  sharply.  This  obser¬ 
vation  could  be  attributed  to  membrane  water  activity  changes 
along  the  channel. 

Fig.  7  shows  the  water  activity  in  the  cathode  flow,  anode  flow 
and  membrane,  respectively.  Near  the  channel  inlet,  the  water 
vapor  in  the  flow  is  sufficient  and  the  membrane  is  well  hydrated, 


Table  2 

Operating  condition  in  base  case 


Parameters 

Values 

Inlet  temperature  of  air 

313  K 

Inlet  temperature  of  hydrogen 

313  K 

Inlet  relative  humidity  of  air 

1.0 

Inlet  relative  humidity  of  hydrogen 

1.0 

Outlet  pressure  of  cathode 

109535  Pa 

Outlet  pressure  of  anode 

109535  Pa 

Current  density 

0.5  A  cm-2 

Excess  coefficient  of  flow  at  cathode 

2.02 

Excess  coefficient  of  flow  at  anode 

1.169 

Fig.  6.  The  distribution  of  current  density  along  channel  in  the  base  case. 


which  increases  the  local  conductive  and  electro-osmotic  drag 
coefficient  of  the  membrane.  As  a  result,  more  hydrogen  ions 
can  pass  through  the  membrane  and  generate  higher  current 
density.  Further  down  the  channel,  the  water  activity  at  the 
anode  flow  drops  quickly  and  the  membrane  becomes  drier  and 
more  resistive  which  decreases  the  current  density.  However, 
since  there  is  some  water  produced  by  the  electrochemical  reac- 


Fig.  7.  The  distribution  of  water  activity  along  channel  in  the  base  case. 
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Fig.  8.  The  distribution  of  hydrogen,  oxygen  and  nitrogen  along  channel  in  the 
base  case. 

tion  continually  at  the  cathode,  the  water  activity  of  the  stream 
increases  gradually  along  the  channel.  Correspondingly,  the 
current  density  increases.  In  summary,  the  water  vapor  fraction 
in  the  flow  has  a  direct  impact  on  fuel  cell  performance.  Fig.  9 
also  displays  that  the  water  activity  in  the  membrane  depends 
on  the  water  activity  in  both  the  cathode  and  anode  streams. 

Fig.  8  shows  how  the  mole  number  of  oxygen  and  hydrogen 
decrease  along  the  channel.  Since  the  excess  coefficient  of  air 
is  higher  than  that  of  hydrogen,  the  molar  fraction  of  oxygen  is 
larger  than  that  of  hydrogen  at  the  channel  exit.  The  reason  for 
choosing  a  large  excess  coefficient  of  air  is  that  the  excess  air 
is  needed  to  take  away  the  extra  water  in  the  cathode  channel. 
The  mole  number  of  nitrogen  does  not  change  because  it  is  not 
involved  in  the  electrochemical  reaction. 

Fig.  9  shows  the  relationship  of  water  activity,  water  con¬ 
tent  and  water  relative  humidity  in  channels  which  have  been 
defined  in  the  model  description  section  above.  In  the  cathode 


Fig.  9.  The  water  activity,  water  content  and  relative  humidity  along  channel  in 
the  base  case. 


channel,  the  relative  humidity  of  stream  is  equal  to  1 .0.  Because 
the  water  is  produced  continually,  the  water  activity  and  water 
content  keep  increasing.  When  the  water  vapor  partial  pressure 
is  greater  than  the  saturation  pressure,  water  vapor  will  condense 
to  liquid  water.  If  the  water  vapor  condensation  rate  is  too  low, 
the  water  vapor  partial  pressure  can  be  greater  than  the  satu¬ 
ration  pressure  in  a  short  period  of  time.  Therefore,  the  water 
activity  can  be  greater  than  1.0.  According  to  their  definitions, 
the  difference  between  the  water  activity  curve  and  water  con¬ 
tent  curve  at  some  point  indicates  that  there  is  liquid  water  in  this 
part  of  the  channel.  In  the  anode  channel,  since  the  flow  is  usu¬ 
ally  unsaturated  and  there  is  no  liquid  water  along  the  channel, 
the  three  curves  are  overlapped.  Fig.  10  shows  the  relationship 
between  water  content  and  water  activity  when  the  anode  inlet 
water  content  is  1.25  and  the  liquid  water  is  injected  into  the 
anode  channels.  In  this  special  case,  although  relative  humidity 
along  the  channels  is  no  more  than  1 .0,  it  still  can  be  seen  that  the 
water  content  curve  and  water  activity  curve  are  not  overlapped. 
The  area  between  water  content  and  water  activity  means  that 
liquid  water  exists  in  the  channels.  This  phenomenon  is  due  to 
the  lower  liquid  water  evaporation  rate.  Meanwhile,  it  is  found 
that  the  water  content  curve  and  water  activity  curve  overlap 
after  60%  of  the  channel  length.  This  behavior  is  attributed  to 
the  fact  that  all  of  the  liquid  water  evaporates  into  water  vapor 
at  this  location.  Therefore,  beyond  this  point,  there  is  no  liquid 
water  in  the  channel. 

Fig.  11  shows  the  temperature  distribution  of  the  cathode 
stream,  the  anode  stream  and  the  stack  along  the  channel.  Fig.  12 
gives  the  detailed  temperature  curves  in  the  vicinity  of  the  inlet. 
At  this  part  of  the  channel,  heat  can  be  transferred  from  the 
stack  to  the  environment  by  convection,  which  leads  to  the  stack 
temperature  being  lower  than  the  stream  temperature.  Further 
down  the  channel,  there  are  several  heat  transfer  processes 
taking  place:  (a)  a  chemical  reaction  occurs  and  reaction  heat  is 
released  to  the  solid  stack;  (b)  water  vapor  condenses  and  latent 
heat  is  released;  (c)  convection  heat  transfer  occurs  as  well  due 
to  the  temperature  difference  between  the  stream  and  stack. 


Fig.  10.  The  water  activity,  water  content  along  channel  for  the  liquid  water 
injection  case. 
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Fig.  11.  The  distribution  of  temperature  along  channel  in  the  base  case. 


Fig.  13.  Distribution  of  pressure  along  channels  in  the  base  case. 


According  to  the  curve,  the  anode  temperature  drops  quickly 
near  the  channel  inlet.  When  it  reaches  the  stack  temperature, 
it  increases  with  stack  temperature.  The  cathode  temperature 
drops  slowly,  at  15%  of  the  channel  length,  it  reaches  the  same 
temperature  as  the  stack.  It  can  be  observed  that  the  large  flowrate 
leads  to  the  flow  temperature  changing  slowly  at  the  cathode. 

Fig.  13  shows  the  distribution  of  pressure  along  the  channel. 
Obviously,  the  pressure  drop  at  the  cathode  channel  is  larger 
than  that  at  the  anode  channel.  This  is  because  the  flowrate  at 
cathode  increases,  which  increases  the  flow  velocity  as  well.  On 
the  other  hand,  the  flowrate  and  velocity  of  the  stream  at  the 
anode  decreases,  consequently,  the  pressure  drop  decreases. 

4.2.  Influence  of  the  inlet  temperatures  of  reactant  gases 

In  this  section,  the  effects  of  inlet  temperatures  of  flow  on  the 
PEM  fuel  cell  performance  will  be  evaluated.  The  inlet  temper¬ 
atures  of  flows  at  both  anode  and  cathode  are  set  to  303,  313, 
323,  333,  and  343  K,  respectively. 


Fig.  12.  The  distribution  of  temperature  at  the  inlet  of  channel  in  the  base  case. 


The  distribution  of  current  density  with  different  inlet 
temperatures  is  shown  in  Fig.  14.  The  distributions  of  current 
density  are  totally  different  for  each  inlet  flow  temperature. 
When  the  inlet  temperatures  of  flow  are  303  and  313  K,  the 
values  of  current  densities  are  highest  at  the  channel  inlet.  They 
then  decrease  quickly  until  they  reach  the  lowest  point  at  around 
the  8%  of  length  down  the  channel.  After  that,  the  current 
density  increases  again.  The  distribution  of  current  density  is 
very  interesting  when  the  inlet  temperature  of  flow  is  323  K. 
The  current  density  increases  slightly  around  the  entrance, 
and  then  begins  to  drop  along  the  channel.  At  approximately 
10%  of  channel  length,  the  current  density  reaches  the  lowest 
value.  After  that,  the  current  density  starts  to  increase  again. 
This  upward  tendency  stops  at  the  position  of  about  60%  of  the 
channel  length.  When  the  inlet  temperatures  of  flow  are  333  and 
343  K,  the  current  densities  increase  near  the  inlet,  and  then  keep 
decreasing  until  the  channel  exit.  This  occurs  primarily  because 
the  current  density  depends  on  the  water  activity  in  the  mem- 


Fig.  14.  A  comparison  of  current  profiles  along  the  channels  with  the  different 
inlet  stream  temperatures. 
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Fig.  15.  A  comparison  of  membrane  water  activity  along  the  channels  with  the 
different  inlet  stream  temperatures. 

brane  and  the  partial  pressure  of  oxygen  in  the  cathode  stream. 
At  a  low  inlet  temperature,  since  the  gas  carries  little  water 
vapor,  water  activity  is  the  primary  factor  that  affects  the  current 
density. 

Fig.  15  shows  the  distribution  of  water  activity  at  the  dif¬ 
ferent  inlet  flow  temperatures.  It  is  found  that  the  local  water 
activities  in  the  membrane  are  less  than  1.0  when  the  inlet  flow 
temperatures  are  303  and  313  K.  At  a  position  around  the  8% 
of  the  length  down  the  channel,  the  water  activity  in  the  mem¬ 
brane  reaches  the  lowest  point.  This  means  the  membrane  is 
very  dry  and  the  speed  of  electrochemical  reaction  is  slow.  Con¬ 
sequently,  only  a  small  number  of  electrons  is  produced,  which 
leads  to  lower  current  density.  With  increasing  membrane  water 
activity  along  the  channel,  more  hydrogen  ions  pass  through  the 
membrane  and  therefore,  the  current  density  increases.  It  is  also 
noticeable  that  the  membrane  activity  increases  quickly  near 
the  channel  exit.  This  can  be  explained  due  to  the  fact  that  the 
stack  loses  heat  to  the  environment,  which  quickly  lowers  the 
stack  temperature  and  flow  temperatures.  According  to  Eq.  (43), 
the  saturation  pressure  will  drop  and  the  water  activities  will 
increase.  For  the  cases  with  higher  inlet  temperature,  such  as, 
333  and  343  K,  the  gases  carry  more  water  vapor  into  the  chan¬ 
nel.  Fig.  15  shows  that  the  water  activity  in  the  membrane  along 
the  whole  channel  is  greater  than  1.0.  According  to  Eq.  (53),  the 
membrane  conductivity  changes  are  small  when  water  activities 
are  large  enough.  Thus,  the  membrane  is  well  hydrated  and  the 
speed  of  electrochemical  reaction  is  fast.  As  a  result,  more  oxy¬ 
gen  is  consumed  and  the  partial  pressure  of  oxygen  decreases 
quickly  (shown  in  Fig.  16).  This  effect  contributes  to  the  drop 
of  the  current  density  along  the  channel.  When  the  inlet  temper¬ 
ature  is  323  K,  the  current  density  depends  on  both  membrane 
water  activity  and  partial  pressure  of  oxygen.  From  the  channel 
entrance  to  about  60%  of  channel  length,  the  membrane  water 
activity  is  less  than  1.0.  The  current  density  changes  with  the 
increasing  membrane  water  activity.  When  the  membrane  water 
activity  is  greater  than  1.0,  the  membrane  conductivity  does  not 
change  much.  However,  the  partial  pressure  of  oxygen  decreases 
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Fig.  16.  A  comparison  of  partial  pressure  of  oxygen  along  the  channels  with  the 
different  inlet  stream  temperatures. 

quickly  and  consequently,  the  current  density  drops  beyond  60% 
of  the  channel  length. 

Fig.  17  shows  that  increasing  the  inlet  temperature  of  flow 
yields  a  higher  cell  potential.  This  behavior  is  attributed  to  the 
fact  that  the  flow  with  high  temperature  introduces  more  water 
into  the  channel  and  membrane  resistance  decreases  due  to  the 
hydration.  It  is  noticed  that  the  polarization  curves  at  the  inlet 
flow  temperatures  of  333  and  343  K  are  overlapped.  This  is 
because  the  membrane  resistance  remains  basically  constant 
when  the  membrane  is  well  hydrated  (water  activity  is  greater 
than  1.0  along  the  whole  channel). 

Fig.  18  shows  the  stack  temperature  distributions.  It  can  be 
seen  that  the  tendency  of  stack  temperature  is  similar  to  that  of 
current  density.  The  larger  the  current  density  is,  the  more  the 
reaction  heat  is  released  and  the  higher  the  stack  temperature. 
Fig.  19  shows  that  the  inlet  temperatures  of  flow  have  a  great 


Fig.  17.  The  effect  of  inlet  stream  temperatures  on  the  performance  of  a  single 
PEM  fuel  cell. 
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Fig.  18.  The  effect  of  inlet  stream  temperatures  on  stack  temperature  of  a  single 
PEM  fuel  cell. 

effect  on  the  pressure  change  along  the  cathode  channels.  When 
the  inlet  flow  temperature  is  high,  the  speed  of  electrochemical 
reaction  increases  and  more  air  is  required.  Thus,  the  flowrate 
in  the  cathode  channels  increase  and  the  pressure  drop  increases 
as  well. 

4.3.  Influence  of  pressure 

Fig.  20  shows  how  the  pressure  affects  the  performance  of  a 
PEM  fuel  cell  under  various  current  densities.  Over  the  entire 
range  of  the  investigated  current  densities,  a  higher  pressure 
leads  to  higher  performance  of  the  fuel  cell.  However,  the  poten¬ 
tial  difference  between  1  and  2  atm  is  greater  than  that  between 
2  and  3  atm.  Furthermore,  this  effect  is  more  obvious  when  the 
current  density  is  high.  This  is  because  the  high-pressure  streams 
can  bring  more  water  into  the  channel  (shown  in  Fig.  21).  As  a 


Fig.  19.  The  effect  of  inlet  stream  temperatures  on  cathode  pressure  of  a  single 
PEM  fuel  cell. 


Fig.  20.  The  effect  of  pressure  on  the  performance  of  a  single  PEM  fuel  cell. 

result,  the  membrane  is  better  hydrated  and  the  speed  of  chemi¬ 
cal  reaction  increases.  Therefore,  the  fuel  cell  can  generate  more 
power  under  the  high  flow  pressure.  From  the  above  analysis,  a 
conclusion  can  be  drawn  that  a  high  inlet  gas  pressure  has  a  posi¬ 
tive  effect  on  system  performance  of  fuel  cell.  However,  whether 
to  use  the  high  pressure  in  a  real  fuel  cell  design  depends  on  the 
trade  off  between  system  improvement  and  the  cost  of  providing 
compressed  gas. 

4.4.  Influence  of  coolant  temperature 

Generally,  the  fuel  cell  system  has  cooling  equipment  to 
remove  waste  heat  and  keep  the  fuel  cell  working  under  optimal 
conditions.  To  simulate  this  kind  of  situation,  it  is  required  to 
investigate  how  the  coolant  temperature  affects  the  fuel  cell  per¬ 
formance.  To  keep  it  simple,  the  coolant  temperature  is  assumed 
to  be  constant.  In  order  to  study  the  effect  on  heat  removal,  the 
polarization  curves  with  different  coolant  temperature,  which 


Fig.  2E  The  effect  of  pressure  on  membrane  water  activity  of  a  single  PEM  fuel 
cell. 
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Fig.  22.  The  effect  of  coolant  temperature  on  the  performance  of  a  single  PEM 
fuel  cell. 

are  293,  298,  and  303  K,  respectively,  are  plotted  in  Fig.  22.  The 
operating  conditions  of  the  fuel  cell  are:  the  flow  inlet  temper¬ 
ature  at  313  K  and  average  current  density  at  0.45  A  cm-2.  As 
can  be  observed,  the  lower  the  cooling  air  temperature  is,  the 
better  fuel  cell  performance  can  be.  The  reason  behind  this  phe¬ 
nomenon  is  quite  simple.  As  the  coolant  temperature  becomes 
lower,  more  heat  can  be  taken  away  from  the  stack,  which 
reduces  the  stack  temperature.  As  the  stack  temperature  goes 
down,  so  do  the  flow  temperatures  due  to  the  intensified  heat 
transfer  between  them  (shown  in  Fig.  23).  As  the  flow  tempera¬ 
ture  becomes  lower,  the  saturation  pressure  of  water  vapor  drops 
which  leads  to  an  increase  in  the  water  activities.  According  to 
Fig.  24,  the  water  activity  in  the  membrane  reaches  the  highest 
when  the  coolant  temperature  is  at  the  lowest  point.  Therefore, 
the  membrane  is  better  hydrated  and  the  speed  of  electrochem¬ 
ical  increases.  As  a  result,  the  performance  of  the  fuel  cell  is 
improved  as  well. 


Fig.  23.  The  effect  of  coolant  temperature  on  stack  temperature  of  a  single  PEM 
fuel  cell. 


Fig.  24.  The  effect  of  coolant  temperature  on  membrane  water  activity  of  a 
single  PEM  fuel  cell. 


4.5.  Influence  of  anode  inlet  humidification 

Water  starvation  in  the  anode  channel  is  one  of  the  prob¬ 
lems  that  fuel  cell  designers  have  to  face.  As  shown  in  the  base 
case  (Fig.  7),  the  anode  and  its  interface  with  the  membrane 
become  less  hydrated  as  the  flow  travels  along  the  channels. 
This  is  because  the  water  vapor  at  the  anode  is  carried  away 
by  hydrogen  ions  and  transported  into  the  cathode.  One  way  to 
solve  this  problem  is  to  inject  liquid  water  into  the  anode  chan¬ 
nel.  As  the  flow  at  the  anode  becomes  unsaturated,  liquid  water 
will  evaporate  to  replenish  the  water  loss.  Therefore,  the  mem¬ 
brane  is  well  hydrated.  Fig.  25  illustrates  variation  of  the  amount 
of  liquid  water  along  the  channel.  It  shows  that  the  liquid  water 
disappears  at  about  25%  of  the  channel  length  when  the  inlet 
water  content  at  the  anode  is  1.1.  If  the  anode  inlet  water  con¬ 
tent  increases  to  be  1.5,  the  liquid  water  will  exist  in  the  whole 
channel.  This  conclusion  can  be  helpful  in  choosing  the  optimal 
anode  inlet  water  content  during  fuel  cell  design. 
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Fig.  25.  A  comparison  of  liquid  water  vary  along  the  channels. 
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4.6.  Constant  stack  temperature  case  analysis 

In  many  previous  studies,  the  stack  temperature  of  the  PEM 
fuel  cell  is  assumed  to  be  constant.  However,  this  temperature 
is  more  likely  to  be  variable  in  real  operation.  To  investigate 
how  the  assumption  of  the  constant  stack  temperature  affects  the 
accuracy  of  the  simulation,  a  constant  stack  temperature  case 
was  simulated  and  the  results  compared  with  those  discussed 
above.  In  this  section,  the  Uniform  Stack  Temperature  can  be 
referred  to  as  “UST”  while  the  Non-Uniform  Stack  Temperature 
is  labeled  as  “NST”. 

Figs.  26  and  13  show  the  temperature  changes  along  the  chan¬ 
nel  for  UST  case  and  NST  case,  respectively.  In  UST  case,  the 
stream  temperature  at  the  anode  is  the  same  as  the  stack  temper¬ 
ature,  since  no  reaction  occurs  and  no  heat  is  generated.  Most 
of  the  reaction  heat  is  taken  away  by  the  stream  in  the  cathode 
channel.  Obviously,  this  case  is  hardly  true  in  real  situations.  The 
curves  in  Fig.  13  are  more  complex,  since  the  reaction  heat  was 
taken  into  account.  Here,  the  anode  stream  temperature  is  also 
close  to  the  stack  temperature  but  changes  along  the  channel.  As 
the  speed  of  local  reactions  change,  the  local  temperature  of  the 
cathode  channel  changes  considerably  as  well.  In  other  words, 
the  temperature  fields  of  the  UST  case  in  both  the  cathode  and 
anode  channels  are  quite  different  to  those  of  NST  case. 

Fig.  27  shows  the  local  current  density  of  the  base  case  along 
the  channel  for  the  UST  case.  As  in  the  case  of  the  correspond¬ 
ing  curve  in  Fig.  6  for  NST  case,  the  current  density  is  highest 
at  the  channel  inlet.  Then,  it  rapidly  drops  to  its  lowest  point. 
After  that,  current  density  increases  slowly  along  the  channel. 
Unlike  the  NST  case,  there  is  no  sudden  increase  of  current  den¬ 
sity  at  the  channel  exit.  This  can  be  explained  by  water  activity 
changes  which  have  been  presented  in  Fig.  28  for  the  UST  case 
and  Fig.  7  for  the  NST  case.  Contrary  to  the  NST  case,  there 
is  no  large  increase  of  water  activity  at  the  end  of  the  channel 
when  the  stack  temperature  is  constant.  Therefore,  the  current 
density  value  does  not  jump  at  the  channel  exit  since  it  is  directly 
influenced  by  the  water  activity.  In  summary,  the  variable  stack 
temperature  can  significantly  affect  the  flow  field  and  thermody- 


Fig.  27.  Current  density  distribution  along  the  channel  in  the  base  case  of  con¬ 
stant  stack  temperature. 

namic  parameters  inside  the  channel.  Its  influence  on  flow  mode 
and  energy  conversion  efficiency  can  hardly  be  ignored.  Also,  in 
the  PEM  model  simulation,  choosing  a  right  boundary  condition 
assumption  is  very  important. 

The  inlet  temperatures  of  flow  have  an  impact  on  the  PEM  fuel 
cell  performance.  Fig.  29  shows  that  a  higher  inlet  temperature 
of  flow  yields  a  higher  cell  potential.  As  discussed  before,  this  is 
mainly  attributed  to  the  fact  that  the  flow  with  high  temperature 
introduces  more  water  to  the  channel  and  decreases  the  mem¬ 
brane  resistance  due  to  the  hydration.  The  distribution  of  current 
density  with  different  inlet  temperatures  is  shown  in  Fig.  30.  It  is 
known  that  the  current  density  depends  on  both  the  water  activity 
in  the  membrane  and  the  partial  pressure  of  oxygen  in  the  cath¬ 
ode  stream.  Figs.  3 1  and  32  show  water  activity  in  the  membrane 
and  the  partial  pressure  of  oxygen  along  the  channel.  At  low  inlet 
temperature,  since  the  gas  carries  little  water,  the  membrane  is 
dry  and  the  speed  of  electrochemical  reaction  is  slow.  As  a  result, 


Fig.  28.  Water  activity  distribution  along  the  channel  in  the  base  case  of  constant 
stack  temperature. 
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Fig.  29.  Polarization  curve  for  different  stack  temperature. 


Fig.  30.  A  comparison  of  current  profiles  along  the  channels  with  the  different 
inlet  stream  temperatures. 


Fig.  31.  A  comparison  of  water  activity  in  membrane  along  the  channels  with 
the  different  inlet  stream  temperatures. 
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Fig.  32.  A  comparison  of  oxygen  partial  pressure  along  the  channels  with  the 
different  inlet  flow  temperatures. 

a  small  amount  of  oxygen  is  consumed  and  the  partial  pressure 
of  oxygen  is  high.  When  the  stream  moves  down  the  channel, 
more  water  is  produced  and  the  membrane  is  hydrated,  which 
increases  the  water  activity  as  well  as  the  speed  of  electrochemi¬ 
cal  reaction.  Therefore,  the  current  density  increases  in  the  later 
part  of  the  channel.  At  high  inlet  temperature,  sufficient  water 
has  been  carried  by  the  gas  at  the  initial  phase  and  the  speed 
of  electrochemical  reaction  is  fast.  As  a  result,  more  oxygen  is 
consumed  and  the  partial  pressure  of  oxygen  decreases  quickly. 
Eventually,  these  factors  lead  to  reduced  electrochemical  reac¬ 
tion  rates  and  decrease  of  the  current  density  throughout  the  rest 
of  the  channel.  Since  the  simulations  provide  different  flow  and 
temperature  fields  inside  the  fuel  cell  channels  for  the  UST  and 
NST,  the  reaction  speed  and  the  amount  of  water  produced  are 
quite  different  at  each  point  in  the  channel.  Figs.  29  and  17  are  the 
polarization  curves  for  the  UST  case  and  NST  case,  respectively. 
In  Fig.  29,  the  difference  between  high  temperature  curves,  such 
as  333  and  343  K,  is  more  apparent  than  for  the  low  temperature 
(313  and  323  K).  On  the  contrary,  the  low  temperature  difference 
is  larger  than  high  temperature  in  Fig.  17.  This  means,  assum¬ 
ing  a  constant  stack  temperature  not  only  affects  the  analysis  of 
flow  conditions,  but  also  changes  the  simulation  result  of  fuel 
cell  performance. 

Since  a  variable  stack  temperature  more  likely  occurs  in  a  real 
fuel  cell,  it  makes  sense  to  replace  the  constant  stack  temperature 
assumption  with  a  variable  stack  temperature.  Actually,  this  was 
shown  in  the  experimental  data.  In  Fig.  6,  the  experimental  result 
has  been  compared  with  both  the  UST  case  and  NST  case.  It  is 
obvious  that  the  NST  case  provides  a  better  simulation  result. 

5.  Conclusions 

In  this  research,  a  model  for  a  single  PEM  fuel  cell  has  been 
developed.  The  simulation  based  on  this  model  can  be  used  to 
analyze  water  transport  across  the  membrane,  the  water  phase 
change  effect,  the  pressure  variation  along  the  channel  and  the 
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energy  balance.  It  can  also  be  used  to  predict  the  characteristics 
of  the  flows  inside  the  channel  and  analyze  the  factors  that  affect 
the  fuel  cell  performance.  Based  on  this  study,  the  following 
conclusions  can  be  drawn: 

1.  The  NST  model  is  more  accurate  than  the  UST  model  when 
predicting  a  single  fuel  cell  performance. 

2.  The  humidification  of  both  anode  and  cathode  sides  is  a  very 
important  factor  affecting  the  performance  of  a  PEM  fuel 
cell. 

3.  Increasing  the  flow  inlet  temperatures  is  an  approach  to  over¬ 
come  the  water  starvation  problem.  However,  if  additional 
equipment  is  added,  the  cost  of  the  fuel  cell  needs  to  be  con¬ 
sidered  as  well. 

4.  Increasing  the  flow  pressure  can  improve  the  fuel  cell  perfor¬ 
mance. 

5.  Proper  liquid  water  injection  at  the  anode  channel  inlet  can 
be  useful  in  fuel  cell  performance  improvement.  An  opti¬ 
mal  amount  of  liquid  water  could  be  determined  by  using 
the  simulations  based  on  the  model  developed  in  the  present 
study. 

6.  Decreasing  the  cooling  temperature  is  helpful  in  improving 
the  fuel  cell  performance. 
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